Primero hay que poder hacer la transformada en 2D. Esta función toma un campo (en forma de vector) y hacer la transformada de Fourier o su inversa. Devuelve una lista con las componentes de Fourier para cada combinación de números de onda en x e y.
fft2d <- function(v, x, y, inverse = FALSE) {
field <- as.matrix(reshape2::dcast(data.table(v, x, y),
x ~ y, value.var = "v"))
field.fft <- fft(field[, -1], inverse = inverse)
dimnames(field.fft) <- list(
k = 1:nrow(field.fft) - 1,
l = 1:ncol(field.fft) - 1)
field.fft <- setDT(melt(field.fft, value.name = "f"))
return(with(field.fft, list(k, l, f)))
}La probamos con una función simple:
test.fun <- function(x, y) {
sin(x*2) + cos(y)
}
test.fun.lap <- function(x, y) {
-4*sin(x*2) - cos(y)
}
test <- setDT(expand.grid(
x = seq(0, 2*pi, length.out = 40),
y = seq(0, 2*pi, length.out = 40)
))
test[, z := test.fun(x, y)]
ggplot(test, aes(x, y)) +
geom_contour(aes(z = z, color = ..level..)) +
scale_color_divergent()
test.fft <- test[, setNames(fft2d(z, x, y),
c("k", "l", "z.hat"))]
ggplot(test.fft, aes(k, l)) +
geom_tile(aes(fill = Mod(z.hat))) +
scale_fill_gradient(low = "white", high = "red")
Sale. Las únicas componetnes con amplitud no nula son k = 4 y l = 1 y sus correspondientes armónicos. La amplitud es muy alta porque no está normalizada, pero no importa por ahora.
Vamos a ver que podemos obtener la función original (a menos de algunas constantes que dan vuelta por ahí) invirtiendo la transformada:
test[, z.inverse := Re(test.fft[, setNames(fft2d(z.hat, k, l, inverse = TRUE),
c("x", "y", "z"))]$z)]
ggplot(test, aes(x, y)) +
geom_contour(aes(z = z)) +
geom_contour(aes(z = z.inverse), color = "red")
Esta función anda bien, así que se puede usar para invertir el laplaciano:
Esta función toma un campo (de nuevo, en forma de vector) y le aplica la inversa del laplaciano a partir de la transformada de fourier.
InverseLaplacian <- function(v, x, y) {
a <- 6371000
field <- data.table(x, y, v)
field[, c("k", "l", "v.hat") := fft2d(v, x, y)]
# field[, lp := l*2*pi/a]
# field[, kp := k*2*pi/(a*cos(y*pi/180))]
field[, v.hat := -v.hat/(k^2 + l^2)]
# field[, v.hat := -v.hat/(kp^2 + lp^2)]
# El modo de k = l = 0 me queda infinito. Como no es demasiado importante (es una constante),
# lo igualo a 0.
field[, v.hat := c(0, v.hat[2:.N])]
field[, c("k", "l", "v.inv") := fft2d(v.hat, k, l, inverse = TRUE)]
return(Re(field$v.inv))
}A partir de la misma función de prueba, voy a calcular su laplaciano y luego invertirlo. La función es muy simple y el laplaciano es el -1 por la función.
test[, z.lap := test.fun.lap(x, y)]
test[, z.inv := InverseLaplacian(z.lap, x, y)]
ggplot(test, aes(x, y)) +
geom_contour(aes(z = z)) +
geom_contour(aes(z = z.inv), color = "red")
¡Éxito!. A menos de una constante aditiva y una multiplicativa, obtengo el campo original.
Vamos a aplicarlo a un campo más realista. La altura geopotencial según NCEP.
ncep <- ReadNetCDF("DATA/NCEP/hgt.mon.mean.nc", vars = "hgt")
setnames(ncep, c("level", "hgt"), c("lev", "gh"))
ncep <- ncep[lev == 300 & date == date[1]]
ncep[, gh := gh - mean(gh), by = lat]Para calcular el laplaciano, voy a tener que derivar:
Derivate <- function(x, y, order = 1, bc = "cyclic") {
# Calcula derivada centrada 1da o 2da
# Entra:
# x: la variable a derivar
# y: la variable que deriva
# order: orden de la derivada (1 o 2)
# bc: condiciones de borde (por ahora, sólo impliementadas cíclicas o nada)
# Sale:
# un vector de la misma longitud que x e y con la derivada
# ¡Asume que la grilla es uniforme!
library(data.table)
N <- length(x)
d <- y[2] - y[1]
if (order == 1) {
dxdy <- (x[c(2:N, 1)] - x[c(N, 1:(N-1))])/(2*d)
} else if (order == 2) {
dxdy <- (x[c(2:N, 1)] + x[c(N, 1:(N-1))] - 2*x)/d^2
}
if (bc != "cyclic") {
dxdy[c(1, N)] <- NA
}
return(dxdy)
}Listo, ahora a partir del laplaciano del geopotencial:
a <- 6731000
ncep[, gh.dxx := Derivate(gh, lon*pi/180, 2)/(a*cos(lat*pi/180))^2, by = lat]
ncep[, gh.dyy := Derivate(gh, lat*pi/180, 2, bc = "none")/a^2, by = lon]
ncep[, gh.lap := gh.dxx + gh.dyy]
ggplot(ncep[abs(lat) < 75], aes(lon, lat)) +
stat_fill_contour(aes(z = gh)) +
geom_contour(aes(color = ..level.., z = gh.lap)) +
map.world.3 +
scale_fill_divergent() +
scale_color_divergent()
Geopotencial en sombreado, laplaciano del geopotencial en contornos.
Puedo intentar obtener el geopotencial original:
ncep[!is.na(gh.lap), gh.inv := InverseLaplacian(gh.lap, lon*pi/180, -lat*pi/180)]
ggplot(ncep[abs(lat) < 80], aes(lon, lat)) +
stat_fill_contour(aes(z = gh)) +
geom_contour(aes(y = lat, z = gh.inv, color = ..level..)) +
map.world.3 +
scale_fill_divergent() +
scale_color_divergent()
Geopotencial en sombreado, laplaciano invertido en contornos
Más o menos. Los centros coinciden un poco, pero la verdad es que es un porquería.
Además, hay un problema adicional. Si no uso los datos globales para calcular invertir el laplaciano, el restulado es distinto:
ncep[!is.na(gh.lap) & lat < 0, gh.inv.SH := InverseLaplacian(gh.lap, lon*pi/180, -lat*pi/180)]
ggplot(ncep[lat < 0 & lat > -80], aes(lon, lat)) +
stat_fill_contour(aes(z = gh.inv)) +
geom_contour(aes(y = lat, z = gh.inv.SH, color = ..level..)) +
scale_fill_divergent() +
scale_color_divergent()
Laplaciano invertido usando todo el globo en sombreado, laplaciano invertido usando sólo el hemisferio suro en contornos.
Probemos con SPEEDY. Como SPEEDY me tira la función corriente y las compontentes de la velocidad, puedo obtener la función corriente a partir de aplicar el laplaciano inverso a la vorticidad y comparar con la solución correcta.
Antes de eso, un diagnóstico de que todo esté bien.
sp <- ReadNetCDF("DATA/SPEEDY/attm.nc", vars = c("psi", "u", "v"))
sp <- sp[lev == 300 & date == date[1]]sp[, v.dx := Derivate(v, lon*pi/180, 1)/(a*cos(lat*pi/180)), by = lat]
sp[, u.dy := Derivate(u, lat*pi/180, 1, bc = "none")/a, by = lon]
sp[, zeta := v.dx - u.dy]
sp[, psi.dxx := Derivate(psi, lon*pi/180, 2)/(a*cos(lat*pi/180))^2, by = lat]
sp[, psi.dyy := Derivate(psi, lat*pi/180, 2, bc = "none")/a^2, by = lon]
sp[, psi.lap := psi.dxx + psi.dyy]
ggplot(sp, aes(lon, lat)) +
stat_fill_contour(aes(z = psi.lap)) +
geom_contour(aes(z = zeta, color = ..level..)) +
map.world.3 +
scale_fill_divergent() +
scale_color_divergent()## Warning: Removed 192 rows containing non-finite values (stat_fill_contour).## Warning: Removed 192 rows containing non-finite values (stat_contour).
El laplaciano de la función corriente en sombreado y la vorticidad relativa en contornos.
Los contornos siguen tan bien al sombreado que prácticamente no se nota la diferencia. Genial.
Hay que notar, igual, que hay una enorme diferencia en la magnitud. En particular, me da que ∇2ψ ∼ 9.8 × 10−7ζ. Esa constante multiplicativa me indica, me parece, que no estoy haciendo bien las constantes que multiplican a las derivadas.
sp[!is.na(zeta), psi.inv := InverseLaplacian(zeta, lon*pi/180, lat*pi/180)]
ggplot(sp, aes(lon, lat)) +
stat_fill_contour(aes(z = psi)) +
geom_contour(aes(z = psi.inv, color = ..level..)) +
map.world.3 +
scale_fill_divergent() +
scale_color_divergent()## Warning: Removed 192 rows containing non-finite values (stat_contour).
Función corriente original en sombreado, función corriente a partir de invertir el laplaciano en contornos.
De nuevo, hay algunas similitudes generales, pero parecería ser una porquería. Aunque veamos…
Calculo U a partir de esta función corriente.
sp[, u.inv := -Derivate(psi.inv, lat*pi/180), by = lon]
sp[, u.psi := -Derivate(psi, lat*pi/180), by = lon]
ggplot(sp[abs(lat) < 85], aes(lon, lat)) +
stat_fill_contour(aes(z = u.psi)) +
geom_contour(aes(z = u.inv, color = ..level..)) +
map.world.3 +
scale_fill_divergent() +
scale_color_divergent()## Warning: Removed 192 rows containing non-finite values (stat_contour).
U derivada a partir de la función corriente real en sombreado, derivada a partir de la función corriente calculada, en contornos.
Epa… no está taaaan mal. Dentro de todo1 y, de nuevo, haciendo la vista gorda a que las magnitudes no tienen nada que verda algo bastante correcto.
Un problema que tiene el método de Fourier es que mi campo no es periódico en dirección “y”. Voy a probar qué pasa si lo hago periódico reflejando el campo.
InverseLaplacian <- function(v, x, y) {
a <- 6371000
field <- data.table(x, y, v)
dy <- unique(y)[2] - unique(y)[1]
field.mirror <- copy(field)
field.mirror[, y := -y + 2*max(y) + abs(dy)]
field <- rbind(field, field.mirror)
field[, c("k", "l", "v.hat") := fft2d(v, x, y)]
# field[, lp := l*2*pi/a]
# field[, kp := k*2*pi/(a*cos(y*pi/180))]
field[, v.hat := -v.hat/(k^2 + l^2)]
# field[, v.hat := -v.hat/(kp^2 + lp^2)]
# El modo de k = l = 0 me queda infinito. Como no es demasiado importante (es una constante),
# lo igualo a 0.
field[, v.hat := c(0, v.hat[2:.N])]
field[, c("k", "l", "v.inv") := fft2d(v.hat, k, l, inverse = TRUE)]
return(Re(field$v.inv))
}sp[!is.na(zeta), psi.inv := InverseLaplacian(zeta, lon*pi/180, lat*pi/180)]## Warning in `[.data.table`(sp, !is.na(zeta), `:=`(psi.inv,
## InverseLaplacian(zeta, : Supplied 8832 items to be assigned to 4416 items
## of column 'psi.inv' (4416 unused)ggplot(sp, aes(lon, lat)) +
stat_fill_contour(aes(z = psi)) +
geom_contour(aes(z = psi.inv, color = ..level..)) +
map.world.3 +
scale_fill_divergent() +
scale_color_divergent()## Warning: Removed 192 rows containing non-finite values (stat_contour).Función corriente original en sombreado, función corriente a partir de invertir el laplaciano en contornos.
Parece mejor.
sp[, u.inv := -Derivate(psi.inv, lat*pi/180), by = lon]
ggplot(sp[abs(lat) < 85], aes(lon, lat)) +
stat_fill_contour(aes(z = u.psi)) +
geom_contour(aes(z = u.inv, color = ..level..)) +
map.world.3 +
scale_fill_divergent() +
scale_color_divergent()## Warning: Removed 192 rows containing non-finite values (stat_contour).
U derivada a partir de la función corriente real en sombreado, derivada a partir de la función corriente calculada, en contornos.